COVID-19 Exploratory

Please feel free to contribute to the Github repo here. You can start by cloning this repository or contacting me at my email to be added as a collaborator. Could certainly use your help!

knitr::opts_chunk$set(
    error = FALSE,
    message = FALSE,
    warning = FALSE
)
library(tidyverse)
library(rvest)
library(magrittr)
library(lubridate)
library(plotly)

Understanding our data sources

John Hopkins CSSE

Daily confirmed cases, deaths and recoveries broken down by country/region and state/province, when state/province is available.

confirmed_ts <- read_csv(url("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv")) %>%
  gather(date, confirmed, -`Province/State`, -`Country/Region`, -Lat, -Long) %>%
  mutate(date = mdy(date))
deaths_ts <- read_csv(url("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Deaths.csv")) %>%
  gather(date, deaths, -`Province/State`, -`Country/Region`, -Lat, -Long) %>%
  mutate(date = mdy(date))
recovered_ts <- read_csv(url("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Recovered.csv")) %>%
  gather(date, recovered, -`Province/State`, -`Country/Region`, -Lat, -Long) %>%
  mutate(date = mdy(date))

joined <- confirmed_ts %>%
  left_join(deaths_ts) %>%
  left_join(recovered_ts) %>%
  gather(stat, value, confirmed, deaths, recovered) %>%
  arrange(`Country/Region`, `Province/State`, date, stat)

joined %>% mutate_if(is.character, as.factor) %>% summary()

          Province/State         Country/Region       Lat              Long              date                   stat      
 Diamond Princess:  330   US            :40755   Min.   :-41.45   Min.   :-157.86   Min.   :2020-01-22   confirmed:25465  
 Grand Princess  :  330   China         : 5445   1st Qu.: 27.61   1st Qu.: -93.06   1st Qu.:2020-02-04   deaths   :25465  
 Adams, IN       :  165   Canada        : 1815   Median : 37.81   Median : -74.30   Median :2020-02-18   recovered:25465  
 Alabama         :  165   Australia     : 1485   Mean   : 31.81   Mean   : -36.13   Mean   :2020-02-18                    
 Alachua, FL     :  165   France        : 1155   3rd Qu.: 42.33   3rd Qu.:  20.94   3rd Qu.:2020-03-03                    
 (Other)         :50985   United Kingdom:  660   Max.   : 71.71   Max.   : 174.89   Max.   :2020-03-16                    
 NA's            :24255   (Other)       :25080                                                                            
     value        
 Min.   :    0.0  
 1st Qu.:    0.0  
 Median :    0.0  
 Mean   :   66.8  
 3rd Qu.:    0.0  
 Max.   :67798.0  
                  

#' Prepares data.frame for plotting by aggregating, weighting, reordering, 
#' and filtering by location.
#'
#' @param geo_level column name to aggregate (sum) by
#' @param dcr_wts weights to order locations; vector of c(deaths, confirmed, recovered)
#' @param show_top_n show top n locations ordered by dcr_wts
#'
#' @return an aggregated, reordered, and filtered by location data.frame
#'
#' @examples
#' plot_prep("Country/Region", c(100, 1, -1), show_top_n = 25)
#'
#' @export
plot_prep <- function(df, geo_level = "Country/Region", dcr_wts = c(1, 0.01, 0),
                      show_top = 1:10) {
  df %>%
    mutate(location = !!sym(geo_level)) %>%
    # aggregate by location
    group_by(location, date, stat) %>%
    summarise(value = sum(value)) %>%
    ungroup() %>%
    # location order
    group_by(location) %>%
    mutate(
      total_wt =
        value[which(date == max(date) & stat == "deaths")] * dcr_wts[1] +
        value[which(date == max(date) & stat == "confirmed")] * dcr_wts[2] +
        value[which(date == max(date) & stat == "recovered")] * dcr_wts[3]
    ) %>%
    ungroup() %>%
    mutate(total_rank = dense_rank(desc(total_wt))) %>%
    # filter top n
    filter(total_rank %in% show_top) %>%
    # order by rank
    mutate(location = fct_reorder(location, total_rank)) 
}

Cross-country comparison

#' Convenience function to plot country comparisons
plot_country_comps <- function(df) {
  df %>% 
    ggplot(aes(date, value, color=location)) +
    geom_line() + facet_wrap(vars(stat), ncol=1, scales = "free_y") +
    theme(legend.position = "right", legend.title = element_blank()) +
    xlab(element_blank())
}

Top countries

joined %>%
  plot_prep("Country/Region") %>%
  plot_country_comps()

excl China

joined %>%
  filter(!`Country/Region` %in%  c("China")) %>%
  plot_prep("Country/Region") %>%
  plot_country_comps()

excl Italy, Iran, South Korea

joined %>%
  filter(!`Country/Region` %in%  c("China", "Italy", "Iran", "Korea, South")) %>%
  plot_prep("Country/Region") %>%
  plot_country_comps()

excl Spain, France

joined %>%
  filter(!`Country/Region` %in%  c("China", "Italy", "Iran", "Korea, South",
                                   "Spain", "France")) %>%
  plot_prep("Country/Region") %>%
  plot_country_comps()

Comparing growth

#' Convenience function to plot growth comparisons
plot_growth_comps <- function(df) {
  df %>%
    ggplot(aes(date, value)) + geom_line() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    facet_wrap(
      vars(location, stat), scales="free_y", ncol = 3,
      labeller=label_wrap_gen(multi_line=FALSE)) +
    theme(strip.background = element_blank(), legend.position = "None") +
    xlab(element_blank())
}

By country

joined %>%
  plot_prep() %>%
  plot_growth_comps()

China

joined %>%
  filter(`Country/Region` == "China") %>%
  plot_prep("Province/State") %>%
  plot_growth_comps()

US

joined %>%
  filter(`Country/Region` == "US") %>%
  plot_prep("Province/State") %>%
  plot_growth_comps()

Canada

joined %>%
  filter(`Country/Region` == "Canada") %>%
  plot_prep("Province/State") %>%
  plot_growth_comps()

Australia

joined %>%
  filter(`Country/Region` == "Australia") %>%
  plot_prep("Province/State") %>%
  plot_growth_comps()

Total deaths since first death

plot_deaths_since_first <- function(df, days_since, init_deaths = 3,
                                    show_top = 1:10) {
  xlab <- paste("Days since deaths =", init_deaths)
  df %>%
    plot_prep("Country/Region", show_top = show_top) %>%
    group_by(location) %>%
    mutate(total_deaths = value[which(date == max(date) &
                                        stat == "deaths")][1]) %>%
    filter(total_deaths > 0) %>%
    mutate(init_death_date = date[which(value >= init_deaths &
                                           stat == "deaths")][1]) %>%
    mutate(
      !!sym(xlab) := date - init_death_date) %>%
    filter(stat == "deaths" & value > 0 &
             between(!!sym(xlab), 0, days_since)) %>%
    ggplot(aes(!!sym(xlab), value, 
               color=location)) + geom_line() +
    theme(strip.background = element_blank(), legend.position = "right",
          legend.title = element_blank()) +
    ylab("Total deaths")
}

Dashboard version

10 days since first death

joined %>%
  plot_deaths_since_first(10)

30 days since first death

joined %>%
  plot_deaths_since_first(30)

60 days since first death

joined %>%
  plot_deaths_since_first(60)

Double days

Double days by country/stat

Top 10 countries

joined %>%
  plot_prep(dcr_wts = c(100, 0.1, 0), show_top = 1:10) %>%
  group_by(location, stat) %>%
  arrange(date) %>%
  mutate(
    `Days to double` = if_else(
      value < 10, as.difftime("", units="days"),
      date - date[sapply(value, function(x) { which.min(abs(x / 2 - value)) })])
  ) %>%
  ggplot(aes(date, `Days to double`, color=stat)) +
  geom_point(alpha = 0.5) +
  geom_smooth() +
  #facet_wrap(vars(location), ncol = 3) + 
  facet_grid(vars(location), vars(stat), scales="free") + 
  theme(legend.position = "bottom", legend.title=element_blank()) +
  xlab(element_blank())

excl China

joined %>%
  filter(`Country/Region` != "China") %>%
  plot_prep(dcr_wts = c(100, 0.1, 0), show_top = 1:10) %>%
  group_by(location, stat) %>%
  arrange(date) %>%
  mutate(
    `Days to double` = if_else(
      value < 10, as.difftime("", units="days"),
      date - date[sapply(value, function(x) { which.min(abs(x / 2 - value)) })])
  ) %>%
  ggplot(aes(date, `Days to double`, color=stat)) +
  geom_point(alpha = 0.5) +
  geom_smooth() +
  #facet_wrap(vars(location), ncol = 3) + 
  facet_grid(vars(location), vars(stat), scales="free") + 
  theme(legend.position = "bottom", legend.title=element_blank()) +
  xlab(element_blank())

Days to double deaths since nth death

h/t: @loeserjohn

Dashboard version

Country/Region level

n_deaths <- 10
min_deaths <- n_deaths * 2
top_n_countries <- 1:20
trunc_to_n <- 2
geo_level = "Country/Region"

# TODO(roboton): Refactor into prep and plot functions
joined %>%
  
  # plot prep at country level
  plot_prep(geo_level, show_top = top_n_countries) %>%
  
  # limit to deaths
  filter(stat == "deaths") %>%
  
  # drop countries with less than min_deaths
  group_by(location) %>%
  filter(value[which.max(date)] >= min_deaths) %>%
  ungroup() %>%
  
  # get the first date with more than or equal to n_deaths
  group_by(location) %>%
  mutate(nth_death_date = date[which(value >= n_deaths) - 1][1]) %>%
  mutate(`Days since nth death` = date - nth_death_date) %>%
  
  # Drop days before nth deaths date
  filter(`Days since nth death` > 0) %>%
  
  # Compute Days to Double
  group_by(location) %>%
  
  mutate(
    # set days to double
    double_idx = sapply(value, FUN=function(x) { max(which(value <= x/2)) }),
    `Days to double` = date - date[double_idx]
  ) %>%
  filter(!is.na(`Days to double`)) %>%
  ## plot
  # truncate to second longest time series
  group_by(location) %>%
  mutate(max_days_since_nth_death = max(`Days since nth death`)) %>%
  ungroup() %>%
  mutate(max_days_rank = dense_rank(desc(max_days_since_nth_death))) %>%
  filter(`Days since nth death` <= max(`Days since nth death`[max_days_rank == trunc_to_n])) %>%
  group_by(location) %>%
  mutate(time_points = n()) %>%
  ungroup() %>%
  filter(time_points > 2) %>%
  ggplot(aes(`Days since nth death`, `Days to double`, color=location)) +
  geom_line(stat = "smooth", method = "auto", se = F) +
  geom_point(alpha = 0.2, size = 1) +
  xlab(paste("days since deaths =", n_deaths)) -> p
  ggplotly(p)

Long view

n_deaths <- 10
min_deaths <- n_deaths * 2
top_n_countries <- 1:15
trunc_to_n <- 1
geo_level = "Country/Region"

# TODO(roboton): Refactor into prep and plot functions
joined %>%
  #filter(!is.na(`Province/State`) & `Province/State` != `Country/Region`) %>%
   
  # plot prep at country level
  plot_prep(geo_level, show_top = top_n_countries) %>%
  
  # limit to deaths
  filter(stat == "deaths") %>%
  
  # drop countries with less than min_deaths
  group_by(location) %>%
  filter(value[which.max(date)] >= min_deaths) %>%
  ungroup() %>%
  
  # get the first date with more than or equal to n_deaths
  group_by(location) %>%
  mutate(nth_death_date = date[which(value >= n_deaths) - 1][1]) %>%
  mutate(`Days since nth death` = date - nth_death_date) %>%
  
  # Drop days before nth deaths date
  filter(`Days since nth death` > 0) %>%
  
  # Compute Days to Double
  group_by(location) %>%
  
  mutate(
    # set days to double
    double_idx = sapply(value, FUN = function(x) { max(which(value <= x/2)) }),
    `Days to double` = date - date[double_idx]
  ) %>%
  filter(!is.na(`Days to double`)) %>%
  ## plot
  # truncate to second longest time series
  group_by(location) %>%
  mutate(max_days_since_nth_death = max(`Days since nth death`)) %>%
  ungroup() %>%
  mutate(max_days_rank = dense_rank(desc(max_days_since_nth_death))) %>%
  filter(`Days since nth death` <= max(`Days since nth death`[max_days_rank == trunc_to_n])) %>%
  group_by(location) %>%
  mutate(time_points = n()) %>%
  ungroup() %>%
  filter(time_points > 2) %>%
  ggplot(aes(`Days since nth death`, `Days to double`,
             color = location)) +
  geom_line(stat = "smooth", method = "auto", se = F) +
  geom_point(alpha = 0.2, size = 1) +
  xlab(paste("days since deaths >=", n_deaths)) -> p
  ggplotly(p) 

Prototype

min_total <- 10
min_stat <- "deaths"
geo_level <- "Country/Region"
max_days_since <- 20
min_total <- 50


joined %>%
  rename(total = value) %>%
  mutate(location = !!sym(geo_level)) %>%
  # aggregate
  group_by(location, stat, date) %>%
  summarise(total = sum(total)) %>%
  # add cfr
  bind_rows(
    (.) %>%
      select(location, stat, date, total) %>%
      spread(stat, total) %>%
      mutate(cfr = deaths/confirmed) %>%
      gather(stat, total, confirmed, deaths, recovered, cfr) %>%
      filter(stat == "cfr")) %>%
  # get max_total and first_date per location/stat
  group_by(location, stat) %>%
  mutate(max_total = max(total),
         first_date = min(date[total >= min_total])) %>%
  group_by(location) %>%
  # drop earlier dates
  filter(date >= first_date[stat == min_stat]) %>%
  # recenter dates
  mutate(days_since = date - first_date[stat == min_stat]) %>%
  ungroup() %>%
  # calc double_days
  group_by(location, stat) %>%
  mutate(
    double_days = date -
      date[sapply(
        total, FUN = function(x) { max(which(total <= x/2)) })]) %>%
  ungroup() %>%
  gather(value_type, value, total, double_days) %>%
  ## plotting
  # truncate days_since
  filter(days_since <= max_days_since) %>%
  # dropping small locations
  group_by(location) %>%
  filter(max_total[stat == min_stat][1] >= min_total) %>%
  ungroup() %>%
  ggplot(aes(days_since, value, color = location)) +
  geom_line() +
  facet_wrap(vars(stat, value_type), scales = "free", ncol=2,
             labeller = label_wrap_gen(multi_line = FALSE))